import numpy as np
import math
#from matplotlib import pyplot
import sys


def Lboundary(x, lx):
    res=x
    while x>=lx:
        x-=lx
    while x<0:
        x+=lx
    return x

def boundary(pos, lx,ly,lz):
    pos[0]=Lboundary(pos[0],lx)
    pos[1]=Lboundary(pos[1],ly)
    pos[2]=Lboundary(pos[2],lz)
    return pos


def anint(x, lx):
    res=x
    if x>0.5*lx:
        res-=lx
    elif x<-0.5*lx:
        res+=lx
    return res

def distancia2(lbox,pos1, pos2):
    dis2=0
    drx=pos1[0]-pos2[0]
    drx=anint(drx, lbox[0])
    dry=pos1[1]-pos2[1]
    dry=anint(dry, lbox[1])
    drz=pos1[2]-pos2[2]
    #drz=anint(drz, lbox[2])
    drx2=drx*drx
    dry2=dry*dry
    drz2=drz*drz
    dis2=drx2+dry2+drz2
    print 'ddd:',drx,dry,drz,math.sqrt(dis2)
    return  [math.sqrt(dis2),drx,dry,drz]

def RShiftedLJ(lbox,pos1,pos2):
    sigma=0.1
    smed=1.0
    r0=smed-sigma
    [dis,drx,dry,drz]=distancia2(lbox,pos1, pos2)
    r_off=dis-r0
    print "r_off: ",r_off
    frac2 = (sigma/r_off)*(sigma/r_off);
    frac6 = frac2*frac2*frac2;
    frep   = 48.0 * frac6*(frac6 - 0.5) / (r_off * dis);
    print 'AAA',frac2,frac6,frep
    return [frep*drx,frep*dry,frep*drz]

def loadConfig(f,npart,lbox):
    pos=[]
    vel=[]
    force=[]
    diameter=[]
    lx=lbox[0]
    ly=lbox[1]
    lz=lbox[2]
    
    for i in range(npart):
        tx = f.readline()
        tx = tx.replace('\n','')
        val = tx.split(' ')
        aux=[]
        for v in val:
            if v!='':
                aux.append(v)
    
        xpos=boundary([float(aux[0]),float(aux[1]),float(aux[2])], lx,ly,lz)
        pos.append(xpos)
        vel.append([float(aux[3]),float(aux[4]),float(aux[5])])
        force.append([float(aux[6]),float(aux[7]),float(aux[8])])
        diameter.append(float(aux[9]))
        
    return [pos,vel,force, diameter]

name='/home/jonk/eclipse/TCLMagneticTest/system_ini.dat'
npart=200
lbox=[8.8858956, 8.8858956, 8.8858956]

f = open(name,'r')
[pos,vel,force, diameter]=loadConfig(f,npart,lbox)
f.close()
iii=0
jjj=19
[dis2,drx2,dry2,drz2]=distancia2(lbox,pos[iii], pos[jjj])
print pos[iii]
print pos[jjj]
print "dist: ",dis2,drx2,dry2,drz2
[fx,fy,fz]=RShiftedLJ(lbox,pos[iii],pos[jjj])
print fx,fy,fz